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Liquid water can become metastable with respect to its vapor in hydrophobic confinement. The 
resulting dewetting transitions are often impeded by large kinetic barriers. According to macroscopic 
theory, such barriers arise from the free energy required to nucleate a critical vapor tube that spans 
the region between two hydrophobic surfaces - tubes with smaller radii collapse, whereas larger ones 
grow to dry the entire confined region. Using extensive molecular simulations of water between two 
nanoscopic hydrophobic surfaces, in conjunction with advanced sampling techniques, here we show 
that for inter-surface separations that thermodynamically favor dewetting, the barrier to dewetting 
does not correspond to the formation of a (classical) critical vapor tube. Instead, it corresponds to an 
abrupt transition from an isolated cavity adjacent to one of the confining surfaces to a gap-spanning 
vapor tube that is already larger than the critical vapor tube anticipated by macroscopic theory. 
Correspondingly, the barrier to dewetting is also smaller than the classical expectation. We show 
that the peculiar nature of water density fluctuations adjacent to extended hydrophobic surfaces - 
namely, the enhanced likelihood of observing low-density fluctuations relative to Gaussian statis¬ 
tics - facilitates this non-classical behavior. By stabilizing isolated cavities relative to vapor tubes, 
enhanced water density fluctuations thus stabilize novel pathways, which circumvent the classical 
barriers and offer diminished resistance to dewetting. Our results thus suggest a key role for fluctu¬ 
ations in speeding up the kinetics of numerous phenomena ranging from Cassie-Wenzel transitions 
on superhydrophobic surfaces, to hydrophobically-driven biomolecular folding and assembly. 


The favorable interactions between two extended hy¬ 
drophobic surfaces drive numerous biomolecular and col¬ 
loidal assemblies m, and have been the subject of 
several theoretical, computational, and experimental in¬ 
quiries [6H22j . Examples include the association of small 
proteins to form multimeric protein complexes, of am- 
phiphlic block copolymers, dendrimers, or proteins to 
form vesicular suprastructures, and of patchy colloidal 
particles into complex crystalline lattices [231127] . When 
two such hydrophobic surfaces approach each other, wa¬ 
ter between them becomes metastable with respect to 
its vapor at a critical separation, d c , that can be quite 
large [8[ |9j 281130) . For nanometer-sized surfaces at am¬ 
bient conditions, d c is proportional to the characteris¬ 
tic size of the hydrophobic object, whereas for micron¬ 
sized and larger surfaces, d c ~ 1 fi m [23 0!]. How¬ 
ever, due to the presence of large kinetic barriers sep¬ 
arating the metastable wet and the stable dry states, 
the system persists in the wet state, and a dewetting 
transition is triggered only at much smaller separations 

(^ 1 nm) [H [221 [2H1 ED]. 

To uncover the mechanism of dewetting, a number of 
theoretical and simulation studies have focused on the 
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thermodynamics as well as the kinetics of dewetting in 
the volume between two parallel hydrophobic surfaces 
that are separated by a fixed distance, d < d c push 
[T6] [HH2J • These studies have highlighted that the bot¬ 
tleneck to dewetting is the formation of a roughly cylin¬ 
drical, critical vapor tube spanning the region between 
the surfaces mi nans], a barrier in the free energetics 
of vapor tube formation as a function of tube radius is 
also supported by macroscopic interfacial thermodynam¬ 
ics, wherein the barrier arises primarily from a compe¬ 
tition between the favorable solid-vapor and unfavorable 
liquid-vapor surface energies (Equation [l] and Figure [l]). 
Thus, the classical mechanism for the dewetting transi¬ 
tion prescribes that a vapor tube, which spans the volume 
between the two surfaces must first be nucleated, and if 
the vapor tube is larger than a certain critical size, it will 
grow until the entire confined volume is dry [9]. 

While it has been recognized that water density fluc¬ 
tuations mm must play a crucial role in nucleating 
vapor tubes, the precise mechanism by which these tubes 
are formed is not clear. To understand how vapor tubes 
are formed and to investigate their role in the dewetting 
process, here we use molecular simulations in conjunc¬ 
tion with enhanced sampling methods mm to char¬ 
acterize the free energetics of water density fluctuations 
in the region between two nanoscopic hydrophobic sur¬ 
faces. Such a characterization of water density fluctua¬ 
tions in bulk water and at interfaces has already provided 
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much insight into the physics of hydrophobic hydration 
and interactions 0H31ISUS3]. In particular, both sim¬ 
ulations and theory have shown that the likelihood of 
observing low density fluctuations adjacent to extended 
hydrophobic surfaces is enhanced relative to Gaussian 
statistics Q31 ED EMM S2]. We show that such en¬ 
hanced water density fluctuations influence the pathways 
to dewetting in hydrophobic confinement by stabilizing 
isolated cavities adjacent to one of the confining surfaces 
with respect to vapor tubes. As the density in the con¬ 
fined region is decreased, the stability of isolated cavities 
relative to vapor tubes also decreases, and at a particu¬ 
lar density, isolated cavities abruptly transition to vapor 
tubes. Surprisingly, for d < d c , that is, separations for 
which dewetting is thermodynamically favorable, we find 
that the nascent vapor tubes formed from the isolated 
cavities are already larger than the corresponding criti¬ 
cal vapor tubes predicted by classical theory. Because the 
newly formed vapor tube is super-critical, it grows spon¬ 
taneously. Importantly, because the formation of this 
super-critical vapor tube involves a non-classical path¬ 
way that circumvents the critical vapor tube altogether, 
the process entails a smaller free energetic cost. Our re¬ 
sults thus point to smaller kinetic barriers to dewetting 
than predicted by macroscopic theory. 



Figure 1 . (a - c) Simulation snapshots of water (shown in 
red/white) in confinement between two square hydrophobic 
surfaces (shown in cyan) of size L — 4 nm that are separated 
by a distance of d — 20 A; configurations highlighting (a) 
the liquid basin, (b) a cylindrical vapor tube of radius, r, 
that spans the confined region, and (c) the vapor basin are 
shown. In the front views, only one of the confining surfaces 
is shown, (d) Macroscopic theory predicts a free energetic 
barrier to vapor tube formation (Equation]!]), suggesting that 
a vapor tube larger than a critical size must be nucleated 
before dewetting can proceed. 


I. MACROSCOPIC THEORY 

According to classical interfacial thermodynamics, the 
free energy for creating a cylindrical vapor tube of radius 
r, which spans the volume between two surfaces sepa¬ 
rated by a distance d, is given by: 

AG t h(A; d) = 7r[r 2 dAP + 2rdy + 2 r 2 7 cos# + 4rA], (1) 

where A P is the difference between the system pressure 
and the saturation pressure, 7 is the liquid-vapor surface 
tension, 6 is the contact angle, and A is the line ten¬ 
sion. For nanoscopic surfaces, the pressure-volume con¬ 
tribution is negligible at ambient conditions |2S102 051, 
whereas the line tension contribution can be impor¬ 
tant [2011351- The term containing cos# is negative for 
hydrophobic surfaces and favors dewetting, whereas the 
term corresponding to formation of the vapor-liquid area 
is unfavorable. The functional form of AGth(j;d) given 
in Equation [l] is illustrated in Figure [l]l for three d-values. 
In each case, a barrier separates the liquid (r = 0) and 
vapor (large r) basins, supporting the notion of dewet¬ 
ting mediated by the nucleation and growth of a vapor 
tube; both the critical vapor tube radius and the barrier 
height increase with increasing d. 

II. DENSITY DEPENDENT FREE ENERGY 
PROFILES FEATURE A KINK 

To investigate how water density fluctuations influence 
the mechanism of dewetting in hydrophobic confinement, 
here we perform molecular dynamics simulations of wa¬ 
ter confined between two, roughly square, hydrophobic 
surfaces of size L = 4 nm, separated by a distance, d, as 
shown in Figure [lji. Water in the confined region is in 
equilibrium with a reservoir of water, which in turn is in 
coexistence with its vapor (A P = 0) [3T1147] , We char¬ 
acterize the statistics of water density fluctuations in the 
confined volume using Indirect Umbrella Sampling (IN¬ 
DUS) HH [32], that is, we estimate the free energy, AG, 
of observing N water molecules in that volume. The free 
energy, AG(7V;d), thus estimated is shown in Figure 
for a range of separations, d, with the free energy of the 
liquid basin, N = Ai iq , being set to zero in each case. 
Over the entire range of separations considered, the free 
energy profile displays distinctive liquid (high N ) and 
vapor (low N ) basins with barriers separating them. In¬ 
terestingly, the free energy profiles also feature a kink, 
that is, an abrupt change in the slope of A G(N; d) is ob¬ 
served at a particular value of N between the liquid and 
vapor basins; we refer to this value of N as Wink- This 
discontinuity in slope is seen more clearly in the deriva¬ 
tives of the free energy, shown in Figure]^. Small errors 
in AG are amplified if simple finite differences are used 
to evaluate the derivatives; we therefore smooth the free 
energy profiles before evaluating the derivatives. Details 
of the smoothing procedure as well as the unsmoothed 
derivatives are shown in the SI. 
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Figure 2. (a) The simulated free energy profiles, 

/ 3AG(N;d ), as a function of the number of water molecules 
between the surfaces, N , display marked kinks (highlighted 
by circles). Here, /? = 1 /feT, with being the Boltzmann 
constant and T being the temperature. The size of the largest 
error bar is also shown for d — 23 A and N — 400. (b) The 
kinks are also apparent in the smoothed derivatives of the free 
energy profiles, which display a sharp decrease in the vicinity 
of ^kink- 

III. KINK SEPARATES THE VAPOR TUBE 
AND ISOLATED CAVITY ENSEMBLES 

To investigate the significance of the kink in the free en¬ 
ergy, we characterize configurations corresponding to N 
on either side of Wink- We do so by building upon the in¬ 
stantaneous interface method of Willard and Chandler Jl] 
to identify iso-surfaces that encompass the dewetted re¬ 
gions, that is, the regions from which water is absent. 
The details of the method are included in the SI. As il¬ 
lustrated in Figure^ for d = 20 A and N = Wink — 12, 
characteristic configurations with N < Wink contain 
a vapor tube, that is, the dewetted region (in purple) 
clearly spans the confined volume between the two sur¬ 
faces (side view, left). Water molecules (not shown for 
clarity) occupy the entire region between the surfaces 
not shown in purple and are also present outside the 
confinement region. In contrast, for configurations with 
N > Wink, isolated cavities are observed adjacent to 



Figure 3. Instantaneous interfaces encompassing dewet¬ 
ted regions (shown in purple) between the hydrophobic sur¬ 
faces (shown in cyan) separated by d = 20 A highlight the 
presence of (a) a vapor tube for N = Wink — 12, and (b) 
an isolated cavity for N =s Wink + 3. Water molecules not 
shown for clarity, (c) Average of the binary vapor tube indi¬ 
cator function, /i tu be, conditioned on the number of waters in 
confinement being V, displays a sharp transition from 1 to 0 
as N is increased. The color scheme is the same as that in 
Figure [5] The value of N corresponding to (/itube) n = 0.5 
(dashed line) is defined as IVtube- (d) IV tu b e is identical to the 
location of the kink in the free energy profiles, Wink, confirm¬ 
ing that the kink demarcates conformations with and with¬ 
out vapor tubes, (e) Conditional average of the isolated cav¬ 
ity indicator function, (/i ca v)iv, shows a sharp increase in the 
vicinity of Wink (the square symbols correspond to IVtube), 
followed by a gradual decrease at larger N- values, and even¬ 
tually vanishing around N = W q - (f) For the larger d- values, 
(/itube)v tube = (/icav)iv tube = 0.5. However, for the smaller d- 
values, (/icav)iv tube < 0.5, suggesting the possibility of direct 
vapor tube nucleation without isolated cavities as intermedi¬ 
ates. 


one surface or the other; however, as seen in Figure [3] 3 
for N = Wink + 3, the cavities do not span the re¬ 
gion between the surfaces to form vapor tubes. Movies 
corresponding to the configurations shown in Figures^ 
and [ 3)3 can be found in the SI. 

The configurations shown in Figures^ an d H> suggest 
that N = Wink marks the boundary between the vapor 
tube and isolated cavity ensembles. To put this notion 
on a quantitative footing, we define indicator functions, 
/itube and /i C av, which are 1 if a given configuration has 
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a vapor tube or an isolated cavity respectively, and 0 
otherwise. The average value of the indicator function 
(htube )n •> subject to the constraint that the number of 
water molecules in confinement is AT, is shown in Fig¬ 
ure for the entire range of AT- values, and for several 
separations, d. Because the configurations were gener¬ 
ated in the presence of biasing potentials, care must be 
exercised in evaluating the averages shown in Figure]^; 
details of the averaging procedure as well as the crite¬ 
ria employed in the definition of the indicator functions 
can be found in the SI. For a given separation, (h tu b e );v 5 
which is the probability of observing a vapor tube con¬ 
ditional on the number of waters in confinement being 
AT, is a sigmoidal function of AT, decreasing sharply from 
1 at low AT to 0 for high AT. We define AT tube to be the 
value of AT at which (h tu be)N undergoes a sharp transi¬ 
tion; in particular, where (A tu be);v crosses 0.5. As shown 
in Figure [3ji, for the entire range of d -values studied here 
(llA<d<25A), AT tube is equal to AT kink , formalizing 
the notion that the kink in AG(AT; d) demarcates con¬ 
figurations that display vapor tubes and those that do 
not. 

Analogous to (At u be);v, the conditional average (h cav )n 
quantifies the fraction of configurations with AT confined 
waters, which feature isolated cavities. For all separa¬ 
tions, (h cav ) 7 v displays a sharp increase, followed by a 
gradual decrease; see Figure [3 Ip. The sharp increase oc¬ 
curs in the vicinity of AT = AT kink as vapor tubes give way 
to isolated cavities, whereas the gradual decrease corre¬ 
sponds to a crossover from isolated cavities to uniform 
configurations. To specify the location of this crossover, 
we define AT cav to be the value of AT where (h cav )N 
crosses 0.5 (with a negative slope). Despite these com¬ 
mon features in the functional form of (h cav ) 7 v, there are 
subtle differences in (h cav )n at small and large separa¬ 
tions, which nevertheless have interesting consequences. 
For the largest separations, a well-defined plateau at 
(h C av )at = 1 separates the sharp increase in (h cav )N 
and its gradual decrease; this plateau demarcates the 
range of AT- values, which reliably feature isolated cavi¬ 
ties. In contrast, the plateau at (h cav )]y = 1 is absent for 
the smaller separations. Instead, the sharp increase in 
(h C av) AT and its gradual decrease overlap in their range 
of AT- values. This overlap suggests that configurations 
corresponding to AT = AT kink are not limited to those 
with vapor tubes or isolated cavities, but could also be 
uniform. Given that (h tube )Ar = 0.5 at AT = AT tube by 
definition, a value of (h cav )N tuhe < 0.5 would correspond 
to a non-zero likelihood of observing uniform configura¬ 
tions with AT = AT tube = AT kink . As shown in Figure [3f, 
that is indeed the case d < 16 A, suggesting that while 
the nucleation of a vapor tube must proceed through the 
formation of isolated cavities for d > 16 A, vapor tubes 
may be nucleated directly from uniform configurations 
for smaller separations. 


IV. FREE ENERGETICS OF VAPOR TUBES 
AND ISOLATED CAVITIES 

Given the abrupt change in the dewetted morphologies 
at AT = AT kink , we expect the functional form of the free 
energetics for AT > AT kink and AT < AT kink to be differ¬ 
ent. Figures [ 3 J 3 and |3jl collectively show that configu¬ 
rations with AT < AT kink feature a vapor tube spanning 
the confined region, consistent with classical arguments. 
While the vapor tube undergoes extensive shape fluctu¬ 
ations, we find its average shape to be roughly cylin¬ 
drical. Coarse-grained density maps of select configura¬ 
tions reflecting the average vapor tube shape are included 
in the SI. To compare the free energetics of the vapor 
tubes obtained from our simulations to macroscopic the¬ 
ory, we first transform the number of waters in the con¬ 
fined region, AT, to an approximate vapor tube radius, 
r, using the simple relation, tt r 2 /L 2 = (AT liq — AT)/AT liq . 
The values of vapor tube radii thus obtained are con¬ 
sistent with the average radii of the vapor tubes ob¬ 
served in our simulations for AT < AT kink , as shown in 
the SI. Having a one-to-one relation between AT and 
r allows us to transform the simulated free energies, 
AG(AT;d), into r-dependent free energies, AG(r; d) in 
the region r k i nk < r < L/ 2; the corresponding free ener¬ 
gies are shown as symbols in Figure^. The lines are fits 
to the AG(r)-data using the macroscopic expectation, 
AGfit(r) = AGth(^) + 2&BTln(l — 2r/L), where the log¬ 
arithmic term corresponds to the translational entropy of 
the vapor tube. AG(r) is fit separately for each d- value 
and yields values of 7 and A that are reasonable. Values 
of 7 are in the range of 12.2 — 15.8 k^T / nm 2 , compara¬ 
ble to the reported value of 14.5 k^T/nvn 2 for the water 
model that we employ [49]. Our fits yield —A /7 in the 
6.5 — 7.5 A range, in accord with a recently reported 
experimental value of A = —30 pN [46], which yields 
—A /7 = 4.2 A. These agreements are remarkable con¬ 
sidering the simplicity of the model that we employ as 
well as the assumptions that we make (cylindrical vapor 
tube shape, constant surface tension independent of va¬ 
por tube curvature, etc.), suggesting that the energetics 
of the vapor tube are well-described by classical macro¬ 
scopic theory. Further details of our fitting procedure, 
the values of the fit parameters for each d, as well as our 
attempts to fit the simulation data to other reasonable 
expressions of G t h(0 can be found in the SI. 

To investigate the free energetics of the isolated cavity 
ensemble, in Figure^, we focus on AG(AT > AT kink ;d). 
In the liquid basin, that is, in the vicinity of AT = AT liq , the 
free energy (symbols) is parabolic (solid lines), indicating 
that the underlying density fluctuations are Gaussian. 
While AG remains harmonic for AT > AT liq , it crosses 
over to being roughly linear (dashed lines) for AT < AT liq . 
Such a crossover from parabolic to linear has also been 
observed adjacent to single extended hydrophobic sur¬ 
faces mi sum], and corresponds to the interfacial water 
undergoing a collective dewetting transition to expel wa¬ 
ter from a nanometer-sized cavity [Ml [38] Eg. Indeed, as 
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shown in Figure^, the location of the crossover agrees 
well with the corresponding A cav -values (squares) dis¬ 
cussed in the previous section. Interestingly, the slopes of 
the linear fat tails are similar for all d-values (in the range 
of ^kink + 20 < TV < 7Vii q — 50); the dashed lines shown 
in Figure (4 Jd are linear fits with a slope of —0.396 k^T 
per water. Additional details of the fitting procedure and 
the values of the parameters obtained are provided in the 
SI. The difference between the values of Ai iq and the x- 
intercepts of the linear fits is also approximately the same 
for all d-values, and is equal to 30 ± 5 waters. Thus, the 
free energy for forming an isolated cavity of a given size 
(as quantified by the number of waters displaced from 
the confined region, Ai iq — TV), is independent of the sep¬ 
aration between the surfaces, that is, the free energetics 
of isolated cavity formation adjacent to one hydropho¬ 
bic surface are largely unaffected by the presence of the 
other confining surface. In contrast, the free energetics of 
vapor tube formation clearly depend on the inter-surface 
separation, d. As a result, the location of the kink, where 
isolated cavities become metastable with respect to vapor 
tubes, also depends on d. 

V. NON-CLASSICAL DEWETTING 

MECHANISM REDUCES BARRIERS 

Figure summarizes our findings for the dewetting 
mechanism presented thus far; in addition to the simu¬ 
lated AG(A; d) for da 15 A, it highlights the metastable 
branches of the vapor tube and isolated cavity ensemble 
free energies, anticipated from the fits shown in Figure [4] 
It is clear that the system minimizes its free energy at 
all times by staying on the branch with the lower free 
energy; at A = Akink, where the two free energy pro¬ 
files intersect, the system jumps from the isolated cav¬ 
ity to the vapor tube ensemble. Importantly, the non- 
classical path leading up to the formation of nascent va¬ 
por tubes (A > Akink) can result in smaller barriers to 
dewetting than anticipated by classical theory, as shown 
in Figure [HJd. For d = 14A, the classical barrier (crit¬ 
ical vapor tube) appears in the metastable segment of 
the free energy profile. The system thus circumvents the 
classical barrier, and instead adopts the path involving 
isolated cavities, which give way to vapor tubes only at 
A = Akink; these nascent vapor tubes are larger than the 
critical vapor tube, so their subsequent growth is down¬ 
hill in energy. Thus, the barrier to dewetting is the free 
energetic cost for forming these nascent, super-critical 
vapor tubes, which is clearly smaller than the classical 
barrier. In Figure^, we illustrate that the nascent va¬ 
por tubes are not super-critical for all separations; for 
d = 23A, the classical barrier appears in the stable seg¬ 
ment of the free energy profile. Thus, while the kink 
in AG(A; d) again marks the formation of a vapor tube 
for d = 23A, the vapor tube formed is smaller than the 
critical vapor tube, and must grow further in a process 
that is uphill in energy, before dewetting can proceed. 



r (nm) 



Figure 4. (a) /3AG(A < Akink; d) is recast as /3AG(r;d), 

the free energy to form a vapor tube of radius, r. The points 
were obtained from the simulated free energy profiles by em¬ 
ploying the relation 7rr 2 /L 2 = (An q — A)/An q in the region 
rkink < t < L/ 2, and the lines are fits to macroscopic the¬ 
ory. (b) The portion of the free energy corresponding to the 
liquid basin (A > Akink) is parabolic at high A (Gaussian 
fluctuations) but linear at low A (fat tails in water number 
distributions). The crossover is gradual and occurs in the 
vicinity of A cav (square symbols), that is, the value of A for 
which (/i C av )n is 0.5 with a negative slope. The linear regions 
have roughly the same slope for all separations. 

As a result, the barrier to dewetting coincides with that 
predicted by macroscopic theory. 

To uncover the separation at which the system transi¬ 
tions from a supercritical to a subcritical nascent vapor 
tube, in Figure [5]i, we plot Akink and A max as a function 
of d. Here, A max corresponds to the value of A between 
the liquid and vapor basins where the AG(A; d) is the 
highest. For small values of d, A max = Akink, indicating 
that the AG(A; d) is the highest at the kink, consistent 
with the formation of supercritical nascent vapor tubes. 
For larger values of d, there is an additional maximum in 
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Figure 5. (a) The simulated f3AG(N ; d) for d = 15A (solid) is shown along with the expected metastable branches of the free 

energies corresponding to the vapor tube (dot-dashed) and isolated cavity (dashed) ensembles. The metastable branches are 
anticipated from the fits shown in Figure [4] The system minimizes its free energy by localizing to the ensemble with the lower 
free energy, (b) For d — 14 A, the nascent vapor tube formed at the kink is larger than the critical vapor tube anticipated 
by macroscopic theory, and therefore grows spontaneously. As a result, the corresponding barrier to dewetting is smaller than 
that predicted by macroscopic theory, (c) For larger separations, here d — 23 A, the newly formed vapor tube is sub-critical, 
and has to grow larger for dewetting to proceed, (d) Comparison of the location of the kink, iVkink, with the location of the 
barrier between the liquid and vapor basins, TV max - For small d, the barrier (point of highest AG) occurs at the kink, so that 
A^ m ax ~ iVkink- In contrast, for larger d- values, the barrier occurs in the vapor tube segment of the simulated free energy profile 
and corresponds to the classical critical vapor tube, so that N < TVkmk- The transition from non-classical to classical behavior 
occurs in the vicinity of the coexistence separation, d c . 


AG at 7V max < A^kink, suggesting that the newly formed 
vapor tubes are smaller than the corresponding critical 
vapor tubes. Interestingly, the separation at which the 
system transitions from non-classical to classical dewet¬ 
ting barriers is close to the separation at which there is 
coexistence between the liquid and vapor, d c . Thus, for 
the separations with a thermodynamically favorable driv¬ 
ing force for dewetting, the mechanism for dewetting is 
manifestly non-classical, corresponding to the formation 
of super-critical vapor tubes from isolated cavities, and 
requiring a smaller free energetic barrier than anticipated 
by macroscopic theory. 

VI. DISCUSSION AND OUTLOOK 

Our results highlight that water density fluctuations 
play a central role in the pathways to dewetting in hy¬ 
drophobic confinement. Enhanced water density fluctu¬ 
ations in the vicinity of hydrophobic surfaces stabilize 
isolated cavities relative to vapor tubes for N > TVkmk; 
the system resides in the classical vapor tube ensemble 
only for TV < TVkmk- While the free energetics of both 
the isolated cavity and vapor tube ensembles are well- 
described by the number of waters in confinement, V, 
this simple order parameter may not be sufficient to de¬ 
scribe the transition from one ensemble to the other. In¬ 
deed, AG(TV; d) represents a projection of a complex free 
energy landscape onto the single parameter, V, and a 
kink in AG (TV; d) strongly suggests that the transition 
from an isolated cavity to a vapor tube involves order pa¬ 
rameter^) that are orthogonal to TV . While beyond the 
scope of the present work, it will be interesting to uncover 
these additional order parameters that define the tran¬ 


sition state ensemble. It is conceivable that additional 
barriers may present themselves in the parameters that 
are orthogonal to TV. 

Dewetting in nanoscopic hydrophobic confinement 
plays an important role in biology; ranging from the as¬ 
sembly of multimeric proteins and the collapse of the 
hydrophobic protein core, to the vapor-lock gating of 
ion channels and the specific binding of ligands to hy¬ 
drophobic grooves on their binding partners. In particu¬ 
lar, recent work has highlighted the importance of includ¬ 
ing the solvent coordinate, V, in describing the kinetics 
of hydrophobically-driven collapse and assembly [5QH53] . 
Our results show that water density fluctuations stabilize 
non-classical pathways, which reduce the barriers along 
the TV-coordinate, and should therefore enhance the ki¬ 
netics of dewetting-mediated biophysical phenomena. 

Dewetting in hydrophobic confinement is also impor¬ 
tant in a host of non-biological phenomena, ranging from 
heterogeneous nucleation of vapor bubbles and contact 
line pinning, to the Cassie-Wenzel transitions on textured 
surfaces [54]. These phenomena involve intricate confine¬ 
ment geometries, which could result in complex pathways 
involving one or more transitions between various dewet- 
ted morphologies, such as the isolated cavities and vapor 
tubes that we observe here. Fluctuation-mediated path¬ 
ways ought to similarly reduce dewetting barriers associ¬ 
ated with these diverse phenomena. 

VII. MATERIALS AND METHODS 

We simulate SPC-E water in confinement between two 
square hydrophobic surfaces of size L = 4 nm, for a range 
of separations, d (See Figure^), ranging from 11 A to 
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25 A, chosen to span the entire range of ^-values with 
both liquid and vapor basins. The surfaces are com¬ 
posed of 1008 atoms each, and are arranged on a hexag¬ 
onal lattice with a spacing of 1.4 A. The surface atoms 
interact with the water oxygens through the Lennard- 
Jones potential with the parameters, a = 3.283 A and 
e = 0.121 kJ/mol; see refs. [201 l2Tj for further details. 
These interactions result in a water droplet contact an¬ 
gle 6 « 120° on the hydrophobic surface, as shown in 
the SI. We have chosen the SPC-E model of water [55] 
since it adequately captures the experimentally known 
features of water such as surface tension, compressibility, 
and vapor-liquid equation of state near ambient condi¬ 
tions, which are important in the study of hydrophobic 
effects mm- Our simulations contain roughly 10,000 to 
15,000 water molecules and were performed in the NVT 
ensemble, thermostatted at T = 300 K using the canon¬ 
ical velocity rescaling thermostat [57] . We employ a pe¬ 
riodic simulation box with the hydrophobic surfaces of 
interest fixed at the center of the box, and a buffering 
liquid-vapor interface nucleated at the top of the box 
with the help of a wall of purely repulsive particles. The 
buffering interface ensures that the system is at the sat¬ 
uration pressure of SPC-E water at 300K Emm; free 


energies obtained with such a construct have been shown 
to be nearly indistinguishable from those obtained in 
the NPT ensemble at a pressure of 1 bar [32] . Short- 
ranged interactions were truncated at 1 nm, whereas 
long-ranged electrostatic interactions were computed us¬ 
ing the particle mesh Ewald method [58]. The bonds in 
water were constrained using SHAKE p59] . To study the 
free energetics of dewetting, we select the cuboid shaped 
(Lx Lxd) observation volume between the hydrophobic 
surfaces, and estimate the free energies, AG(N ; d), using 
the indirect umbrella sampling (INDUS) method |31,32]. 
Each biased simulation was run for 6 ns and the first 1 ns 
was discarded for equilibration. 
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I. ESTIMATING “SMOOTHED” DERIVATIVES 
OF THE FREE ENERGY 


II. OBTAINING INSTANTANEOUS 
INTERFACES THAT ENVELOP DEWETTED 
REGIONS 


To demonstrate the existence of kinks in A G(N; d ) 
more clearly, we plot smoothed derivatives of A G(N; d ) 
in Figure 2b of the main text. Directly computing deriva¬ 
tives of AG (TV) using finite difference, for example, using 
dAG/dN = A G(N + 1) — AG(TV), results in noisy esti¬ 
mates due to the numerical error in A G(N) (Figure [Si]). 
To minimize such effects, we first smooth A G(N) using 
a rectangular window function, such that 


AG sm0 othed (TV) — 


1 

AN 


N+AN/2-1 

X AG ( 7V ) 

N—AN/2 


(S2) 


is the smoothed free energy. Smoothed derivatives are 
then estimated from finite differences of such smoothed 
free energies, 


^k(TV) » [AG smoothed (N + 1) - AG smoothed (iV)]. 

(S3) 

Derivatives of both the unsmoothed and the smoothed 
(using a width of AN = 10) free energies for select values 
of d are compared in Figure [STJ 


An analogous smoothing of the first derivative was also 
performed in order to obtain the second derivative of the 
free energy, which in turn enables us to robustly estimate 
TVkink as the location of the minimum in d 2 f3AG/dN 2 . 
The error in TVkink was estimated using block averaging; 
the trajectories were divided into five blocks of 1 ns each, 
and TVkink was obtained for each of the five blocks. As an 
example, the second derivatives for five such blocks are 

o 

shown in Figure S2 for d = 14 A . 
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this work 
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Here, we closely follow and build upon the approach for 
estimating instantaneous interfaces, originally developed 
by Willard and Chandler [I]. In addition to consider¬ 
ing the coarse-grained water density, we also include a 
coarse-grained density arising from the plate atoms into 
an overall, normalized coarse-grained density at (x,y, z) 
as follows: 


~ / Pwater (*G IT %•) H) Pplate(*G Vi H) 

Ptotai(ah2/,2;R) = ■ 


Umlk 
r water 


Opiate 


(S4) 

Here, R represents the positions of all the heavy atoms 
in the system in a given configuration, the coarse-grained 
water density is normalized by the corresponding bulk 
density, and the coarse-grained plate density by the max¬ 
imum coarse-grained plate density that occurs at the cen¬ 
ter of the plate. Defined this way, p to tai is approximately 
equal to unity everywhere in the liquid state and near 
unity at the center of the plates and in the interfacial 
region. Configurations containing dewetted regions (cav¬ 
ities) will have significantly smaller values of ptotai that 
approach zero in the vicinity of the cavity. Therefore, 
the ptotai = 0.5 iso-density surface serves as a convenient 
definition of the instantaneous interface, allowing us to 
readily visualize the position, size, and shape of the cav¬ 
ity; we use the Marching Cube algorithm [2j to identify 
the instantaneous interface. The coarse-grained density 
fields of the individual species are estimated using 


_ N a 

Pa(x,y,z; R) = ~ x )^Vi ~ v)H z i ~ z ), (S5) 

i=1 


where a represents either the water oxygen atoms or the 
plate atoms, N a is the number of atoms of type a , and 
(xi,yi,Zi) correspond to the coordinates of atom i. 

For each configuration obtained from our simulations, 
we set up a three-dimensional grid to compute the coarse¬ 
grained density field with 0.1 nm spacing in each dimen¬ 
sion. The coarse-graining function cf(x) is chosen to be 
a Gaussian with a width of 0.24 nm, which was trun¬ 
cated at 0.7 nm, shifted down to make it continuous, 







and normalized. The particular characteristics of the in¬ 
stantaneous interfaces thus computed, as well as those of 
the dewetted regions, such as their exact shapes and vol¬ 
umes, will depend on the choices made in Equations |S4| 
and[S5j as well as the parameters chosen. A discussion of 
how these choices affect the instantaneous interface cal¬ 
culation and which choices are judicious is beyond the 
scope of this work, and will be the subject of a separate 
publication. However, it is important to note that the 
qualitative insights that we obtain in this work are not 
sensitive to the particular choices that we make here. 


III. MOVIES 


The dynamical nature of the vapor tubes described in 
the main text can be observed in the movie shown in 
Figure [S3j There, we show a top and a side view of a 
plate-spanning vapor tube from a simulation with a bi¬ 
asing potential, k(N — TV*) 2 /2 ; k = 0.12 kJ/mol and 
V* = 650. The average value of TV in the presence of the 
biasing potential was ^ N^j = 653. The plate atoms are 

shown as cyan spheres, water molecules are not shown 
for clarity, and the purple mesh corresponds to the in¬ 
stantaneous interface enveloping the vapor region. In 
essence, the purple mesh defines regions devoid of water 
molecules, such that the space outside the mesh is high 
in water density. 

We show similar movies for an isolated cavity at the 
surface(s) of the hydrophobic plate(s) in Figure |S4| The 
coloring scheme is the same as in Figure [S3j and the 
configurations are taken from a biased simulation with 
TV* = 660 and k = 0.12 kJ/mol; the corresponding 

^TV^ = 668. There, an isolated cavity can be observed 

fluctuating in shape and size, and even moving from one 
plate to the other. The mechanism for bubble migration 
from one plate to another is an open question relevant to 
dewetting transitions and deserves further investigation. 


IV. DEFINITION OF THE TUBE AND CAVITY 
INDICATOR FUNCTIONS 


The tube indicator function h tu be is determined by ex¬ 
amining the coarse-grained density field, p to tai> between 
the two plates. We define the x-coordinate to be perpen¬ 
dicular to the plates, such that the two plates are located 
at x\ and £ 2 , respectively, with x\ < x 2 - We additionally 
define a buffer region b = 0.4 nm from the center of each 
plate to avoid the region where the coarse-grained density 
originating from the plate atoms is larger than or close to 
0.5. If at some location, (?/*, z*), the total coarse-grained 
density is below 0.5 at all x -values between x\ + b and 


X 2 — b , we assign htube = 1? i.e. 


fttube = lif 3(2/*,**) 

: {ptotaifa, V*i z *) < °- 5 l V x e Oi + 6, x 2 - b) 
= 0 otherwise. (S6) 

The cavity indicator function h cav is similarly deter¬ 
mined from ptotai- However, h cav is equal to unity if the 
coarse-grained density is less than a half at some location 
between the plates, but not for all x-values between the 
plates. In other words, 

^ cav = 1 if 3 (^*, 2 /*,^*) 

:{/W^V*)<0.5}, h tube = 0 

= 0 otherwise. (S7) 


V. PROCEDURE FOR CALCULATING 
UNBIASED ENSEMBLE AVERAGES 


All umbrella sampling simulations performed in this 
work employ a biasing potential. Therefore, when per¬ 
forming ensemble averages of any observable, care must 
be exercised in accounting for the bias introduced by the 
potential. This is done using the WHAM/MBAR for¬ 
malism |3H5], such that the unbiased ensemble average 
of any observable, A(R), which can be expressed as a 
function of the configuration vector, R, is given by 


K Nj 

( A m=c- i j2E 


3= 1 n=1 


_ 


(S8) 


where 


K Nj 

c = --—=—. (S9) 

Here, Fk is the free energy of the kth biased ensemble, K 
is the number of biasing potentials or windows used, TV& 
is the number of samples in window fc, T4 is the biasing 
potential for window fc, and Ry n refers to configuration 
n in window j. 

Equation [S8| is used here to calculate the ensemble av¬ 
erage of h tu be 5 conditioned on the number of waters in 
confinement being TV, according to 


. , _ . \^tube(R)^AT,W(R) ) 

(^t U be(R))jy = ^—j -(S10) 

\^/V,AT(R) / 

where S N ^ is the Kronecker delta function. The en¬ 
semble average of h cav as a function of TV is calculated in 
an analogous manner by replacing h tu b e (R) with h cav (R) 
in Equation |S10[ 







VI. AVERAGE SHAPE OF THE VAPOR TUBE VIII. DETAILS OF VAPOR BASIN FITS 


The average shape of a vapor tube is identified by 
first taking the average of the coarse grained density, 
Ptot(z,2/,z) = (/Ototal over 5000 frames. Fig- 

o 

displays the averaged density of the d = 20 A 


S5 


ure 

with an average of 569 waters (top) and 474 waters (bot¬ 
tom) between the plates, obtained from biased simula¬ 
tions with TV* = 570 and k = 0.03 kJ/mol and TV* = 480 
and k = 0.03 kJ/mol, respectively. From left to right 
in each row presents a front view, side view and three- 
dimensional view of the averaged shape of the vapor tube. 
The three-dimensional rendering of the vapor tube is ob¬ 
tained as the iso-density surface p t ot = 0.5, and this 
surface provides an accurate description of the average 
shape of the vapor tube. Additionally, the radii obtained 
directly from these isodensity surfaces, 1.03 nm and 
1.30 nm, respectively for the top and bottom panels of 
Figure [S5j are in reasonably good agreement with those 
obtained from the simple re latio n between r and TV used 
in the main text (Equation |S13| ), 1.22 nm and 1.45 nm, 
respectively. Note that the vapor tubes are cylindrical to 
a good approximation, in accord with macroscopic the¬ 
ory. However, further corrections to such theories could 
be obtained by taking into account finer details of the 
“hour-glass” shape of the vapor tubes. Because the pre¬ 
cise shape of the vapor tubes depends on the details of 
the instantaneous interface calculation and the parame¬ 
ters employed, we do not attempt such corrections here. 


VII. CONTACT ANGLE DETERMINATION 


The contact angle of the surface was determined by 
performing a 2 ns simulation of a cylindrical droplet con¬ 
taining 4142 water molecules. This cylindrical droplet 
was divided into five slabs, each 1 nm in width, and the 
density as a function of the radius y and the height z 
of each slab was computed. The droplet profile is then 
defined as the point where the density of the droplet is 
equal to half that of the bulk density. 

The contact angle was then determined by fitting this 
droplet profile to a circle for z > 0.7 nm, as shown in 
Figure [S6| where the fit function is given by 

y = \J&- ( z — a ) 2 - (Sii) 


A. Fitting used in the main text 

In order to fit the simulated free energies to macro¬ 
scopic theory, we consider the formation of a vapor tube 
of radius r, which depends on the number of waters, TV, 
between the hydrophobic plates as 


r 


( N lUl -N \ V 

V Miq ) 7T ' 


(S13) 


Note that Equation mu is an approximate expression; 
however, this simple approximation captures the salient 
features of vapor tube formation, as detailed below and 
in the main text. The free energy as a function of r as 
predicted by macroscopic theory is 


/3AG(r) = 271/3^ 


'l n I 7 2A 

r cos 6 + r d e ff H- 

7 


2 r 

—2 ln(l——), 
(S14) 

where 7 is the liquid-vapor surface tension, 6 is the con¬ 
tact angle (determined as described above), and A is the 
line tension. The first term in Equation |S14| is the free 
energy of vapor tube formation, AGth> described in the 
main text, and the last term accounts for the transla¬ 
tional entropy of the vapor tube. The effective distance 
between the plates, d e ff, is obtained by subtracting a con¬ 
stant offset, £, from d, that is d e q = d—£. The x-intercept 
of a linear fit of the simulated d-dependence of TVi iq is 
equal to £, as shown in Figure [S7j so that a plot of TVi iq 
vs d e ff passes through the origin. Here we find a value 
of the offset to be £ = 0.4964 nm. The fits shown in 
Figure 4a of the main text were obtained by fitting the 
simulated free energies for TV < TVkink and r < Lj 2 to 
the parameters, /J7 and A/7, using Equation |S14| the fit 
parameters obtained are listed in Table [TJ 


B. Curvature Corrections 

We also explored a number of other possible fits, the 
first of which includes curvature corrections to the sur¬ 
face tension using the Tolman length, d, such that the 
curvature corrected surface tension is 7(7*) = 7(1 — S/r). 
The fit equation then becomes 


Using this functional form, the contact angle, d, can be 
obtained from 

dy (z — cl) 1 

dz z=0 7 „ m ~~ [6 2 — (z — a) 2 ] 1/2 " tan(7r - 9 )' 

(S12) 

The contact angle was determined independently for each 
of the five slabs and averaged to yield a contact angle of 
0 = 120.2° ± 0.5°, or cos 0 = -0.503 ± 0.008. 


f3AG(r) = 27T/37 


2A 


7 


2i " 11 -? 


r cos 0 + r ( d e ff H- S cos 6 ) — Sd e q 

(S15) 


However, by referencing the simulated free energies for a 
given d-value to the free energy in the liquid basin, that 
is AG(TVii q ) = 0, we force AG(r) to be 0 at r = 0. To be 
consistent with this convention, we neglect the constant 
















TABLE I. Parameters obtained from fitting the vapor tube 
portion of the free energy to Equation |S14| 


d( A) 

/?7 (nm' 

11 

15.1 

12 

15.4 

13 

15.4 

14 

15.4 

15 

15.7 

16 

15.0 

17 

15.1 

18 

14.6 

19 

14.2 

20 

13.7 

21 

13.1 

22 

13.3 

23 

12.9 

24 

12.5 

25 

11.9 


2 ) A /7 (nm) 
-0.150 

-0.167 

-0.177 

-0.180 

-0.180 

-0.207 

-0.208 

-0.223 

-0.230 

-0.230 

-0.243 

-0.240 

-0.234 

- 0.221 

-0.216 


term, in Equation |S15[ and instead fit to 


13AG(r ) = 27t/?7 


r z cos 0 + r ( d e ff H—— — 5 cos 0 


7 




(S16) 


From Equation |S16| we see that the curvature correc¬ 
tion S acts in a manner analogous to the line tension A, 
i.e. both are coefficients in the term linear in r. There¬ 
fore, we can simply relate the fit parameters in Equa¬ 
tion |S16| to those in Equation |S14[ according to 


X X' 5 
- = — + - cos <9, 

7 7 ' 2 


(S17) 


where the prim ed quantities indicate those obtained from 
Equation S14 The value of A/7 obtained from Equa¬ 
tion [ST6] can then be readily predicted with knowledge of 
the Tolman length. Recent estimates of this length for 
SPC/E water at 300 K yield S ~ —0.1 nm [6], 7j. There¬ 
fore, the value of A/7 changes by ten to sixteen percent 
through the inclusion of the Tolman length in our fitting 
procedure. 


C. Omission of Line Tension 

To ascertain the importance of line tension, we attempt 
to fit the simulated free energies without including the 


term containing A in our expression for the macroscopic 
theory, and instead using S as a fit parameter according 
to 


/ 3AG(r ) = 2tt/3 7 [r 2 cos 6 + r ( d e ff — S cos 0)] 




(S18) 


As one may anticipate from the above discussion, the 
data is fit equally well by Equation |S18| However, 


this procedure yields unphysical estimates of the Tolman 
length, with S ~ — 1 nm. Therefore, we can conclude that 
line tension is necessary to provide a physically accurate 
description of the vapor tube formation and growth that 
facilitates capillary evaporation between nanoscopic hy¬ 
drophobic plates. 


IX. LIQUID BASIN FITS 


In order to characterize the free energetics of isolated 
cavities, we fit the simulated free energies in the range, 
TVkink + 20 < N < TViiq — 50, to straight lines. We first 

o 

fit the d = 25 A free energy, and use the corresponding 
slope for all d-values. Intercepts are then obtained by fit¬ 


ting the data for each d-value separately. For d < 16 A , 
data in the + 20 < N < Nn q — 50-range are insuffi¬ 
cient to be fit reliably, so the fitting was not performed. 
The regions fit for each d-value are shown in Figure |S8| 
In order to estimate errors, we divided the data into five 
blocks of equal length, and fit each block separately (sim¬ 
ilar block averaging analysis was performed for all quan¬ 
tities). Therefore, for each d, we show five data sets and 
five fits, although they are nearly indistinguishable. 

From these fits, we can examine the behavior of 
TViiq(d) — JVint(d), where 7V int (d) is the x-intercept of the 


linear fit to the fat tails (for d > 16 A ). Both Ny iq (d) and 


JVint(d) are roughly linear in d, see Figur es |S7| and S8 


re¬ 


spectively. However, as shown in Figure |S9[ the quantity 
TViiq — 7Vj nt , which represents the distance in N from the 
liquid basin one must move to observe non-Gaussian tails 
in A G(N] d), or equivalently, to form isolated cavities, is 
largely insensitive to d. 


X. BARRIER LOCATION AND HEIGHT: 
SIMULATIONS VS MACROSCOPIC THEORY 


In this section, we focus on the location as well as the 
height of the barrier, as predicted by macroscopic theory, 
and compare them with the corresponding simulated val¬ 
ues. Macroscopic theory predicts a linear dependence of 
the location of the barrier, r*, on d, 


r*(d) = — 




A 


(S19) 


2 cos 6 7 cos 0 

where we have neglected the logarithmic term in Equa¬ 
tion S14 This behavior is captured by the data shown in 
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Figure [SlOl where the revalues are obtained from the fits 
of the vapor tubes free energetics to macroscopic theory. 
The slope of the fitted line corresponds to cos 0 = —0.547, 
which is slightly larger in magnitude than that obtained 
from direct simulations. The ratio A/7 = —0.16 nm ob¬ 
tained from the ^/-intercept of the linear fit lies in the 
same range as those predicted from fitting the free ener¬ 
gies to Equation S14 Macroscopic theory also predicts 
AG(r*; d) to be quadratic in d\ this behavior is similarly 
captured as illustrated by fitting the data to a parabola 
(solid line). 

In addition to the position and height of the maximum 
of the fits, we also include the position and the height of 
the free energy barrier obtained from the simulated free 
energy profiles, r max (d) and AG(r max ;d), respectively. 
For large d- values, we expect r max ~ r*, however, this is 


not true at smaller plate separations because the critical 
vapor tube is in the metastable branch of the vapor tube 
free energy; for these d-values, the barrier corresponds to 
the kink in the free energy. Indeed, r* 7^ r max at these 
separations, and the dependence of r max on d can not be 
described by a straight line over the entire range of d. 
Similar behavior is observed in AG(r max ;d), albeit to a 
lesser extent. 


XI. DATA FOR ALL SYSTEMS STUDIED 


For completeness, we include the free energies, their 
derivatives, (h tu b e ) N , an d (^cav)jv> as we ^ as fits to the 
vapor and liquid basins, for all d -values studied in Fig¬ 
ures Sll[ |S12[ |S13[ |S14[ and S15 respectively. 


[1] Willard, A. P.; Chandler, D. Instantaneous Liquid Inter¬ 
faces , J. Phys. Chem. B, 114, 1954-1958 (2010). 

[2] Lorensen, W. E.; Cline, H. E. Marching cubes: A high 
resolution 3D surface construction algorithm , Computer 
Graphics, 21, 163-169 (1987). 

[3] Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte 
Carlo data analysis , Phys. Rev. Lett., 63, 1195-1198 
(1989). 

[4] Tan, Z.; Gallichio, E.; Lapelosa, M.; Levy, R. M. The¬ 
ory of binless multi-state free energy estimation with ap¬ 
plications to protein-ligand binding , J. Chem. Phys., 136, 
144102 (2012). 


[5] Shirts, M. R.; Chodera, J. D. Statistically optimal analy¬ 
sis of samples from multiple equilibrium states , J. Chem. 
Phys., 129, 124105 (2008). 

[6] Sedlmeier, F.; Netz, R. R. The spontaneous curvature 
of the water-hydrophobe interface , J. Chem. Phys., 137, 
135102 (2012). 

[7] Vaikuntanathan, S.; Geissler, P. L. Putting water on a 
lattice: The importance of long wavelength density fluctu¬ 
ations in theories of hydrophobic and interfacial phenom¬ 
ena. , Phys. Rev. Lett., 112, 020603 (2014). 

[8] Sharma, S.; Debenedetti, P. G. Evaporation rate of water 
in hydrophobic confinement Proc. Natl. Acad. Sci. USA, 
109, 4365-70 (2012). 








14 



Figure SI. Derivatives of the unsmoothed (points) and 
smoothed (lines) free energies for d — 14 A (blue) and d = 
23 A (red). 



400 425 450 
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Figure S2. Second derivative of the free energy for d = 

o 

14 A . Block averaging was used to obtain iVkink and its 
associated error bars; the derivatives of each of the the five 
blocks are shown here. iVkink is the location of the minimum 
in d 2 /3AG/dN 2 . 
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Figure S3. Movies (click to view) illustrating a plate- 
spanning vapor tube in a biased simulation with N* = 650 
and Hi — 0.12 kJ/mol, with — 653. Plate atoms are 

shown as spheres, while the purple mesh corresponds to the 
instantaneous interface enveloping the vapor region (water 
molecules have been omitted for clarity). The movie is shown 
from directions (top) perpendicular and (bottom) parallel to 
the plane of the plates. The duration of the movie corre¬ 
sponds to 250 ps of simulation time. Adobe Reader may be 
needed to view. Representative snapshots are shown here; 
movie files are available upon request and will appear in the 
final published version. 
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Figure S4. Movies (click to view) illustrating vapor bubbles 
formed at the surface of a hydrophobic plate in a biased simu¬ 
lation with N* = 660 and k = 0.12 kJ/mol, with — 668. 
Plate atoms are shown as spheres, while the purple mesh cor¬ 
responds to the instantaneous interface enveloping the vapor 
region (water molecules have been omitted for clarity). The 
movie is shown from directions (top) perpendicular and (bot¬ 
tom) parallel to the plane of the plates. The duration of 
the movie corresponds to 250 ps of simulation time. Adobe 
Reader may be needed to view. Representative snapshots are 
shown here; movie files are available upon request and will 
appear in the final published version. 
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Figure S5. The average coarse-grained density, ptot(x, y , z), 

o 

for d = 20 A , with an average of 569 waters (top) and 474 wa¬ 
ters (bottom) between the plates. The top panel was obtained 
from a biased simulation with N* = 570 and k = 0.03 kJ/mol, 
and the bottom panel was obtained from a biased simulation 
with TV* = 480 and k = 0.03 kJ/mol The rightmost panels 
depict three-dimensional renderings of the vapor tube shape 
at the iso-density surface p to t = 0.5 

















18 




0.5 1 1.5 2 2.5 3 3.5 

y (nm) 


Figure S6. (top) Snapshot of the cylindrical droplet used 
to calculate the contact angle, (bo ttom ) Fit of the average 
droplet profile (points) to Equation |S11| (line). 
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Figure S7. The number of water molecules between the 
plates in the liquid basin, iVii q , varies linearly with the dis¬ 
tance between the plates, d, as measured between the centers 
of the atoms of the two plates. The effective distance between 
the plates, d e ff, is obtained subtracting from d the x- intercept, 
£, of a linear fit of iVn q (d); thus, d e ff = d — £.. 
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Figure S8. (top) Portions of the free energies that were 

o 

linearly fit for d > 16 A . Simulation data, split into five 
blocks for error estimation, is shown as colored data points. 
Linear fits are shown as black dashed lines, (bottom) The 
intercept, ATit, obtained from the linear fits, is plotted as a 
function of d. 
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Figure S9. The number of waters needed to be removed 
from the liquid basin to encroach on the region of the fat 

o 

tail, as measured by iViiq — Mnt for d > 16 A . Error bars 
correspond to one standard deviation. 








22 




15 20 25 


d (A) 


Figure S10. (top) Barrier position r* and (bottom) barrier 
height /3AG(r*) of the fitted vapor tube free energies, are well- 
fit by linear and parabolic functions of the inter-plate sepa¬ 
ration, d, respectively, in agreement with macroscopic theory. 
In contrast, the simulated barrier position r max and barrier 
height /3AG(r max ) deviate from those functional forms. 
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Figure Sll. The dependence of (top) the free energies 
and (bottom) their smoothed derivatives on N for all d -values 
studied. Arrows point in the direction of increasing d , from 

o o o 

d — 11 A to d = 25 A in increments of 1 A . 
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Figure S12. (ht u be) N f° r values studied. Arrows point 
in the direction of increasing d, from d = 11 A to d = 25 A 

o 

in increments of 1 A . 



Figure S13. (h ca v) N for all d-values studied. Arrows point 

o o 

in the direction of increasing d, from d — 11 A to d = 25 A 
in increments of 1 A . 
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Figure S14. Free energies in the vapor basin (N < iVkink) 
and the corresponding fits to Equation |S14| Simulation data 
is shown as points and the fits are shown as solid lines. The 

o 

arrow indicates the direction of increasing d, from d = 11 A 
to d = 25 A in increments of 1 A . 



Figure S15. Free energies in the liquid basin (.N > iVkmk) 
and the corresponding parabolic (solid) and linear (dashed) 
fits to the right side of the minimum and to the fat tail regions, 
respectively. 



















